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ABSTRACT 

We have investigated planetary accretion from planetesimals in terrestrial planet 
regions inside the ice line around M dwarf stars through iV-body simulations including 
tidal interactions with disk gas. Because of low luminosity of M dwarfs, habitable zones 
(HZs) are located in inner regions (~ O.lAU). In the close-in HZ, type-I migration and 
the orbital decay induced by eccentricity damping are efficient according to the high disk 
gas density in the small orbital radii. Since the orbital decay is terminated around the 
disk inner edge and the disk edge is close to the HZ, the protoplanets accumulated near 
the disk edge affect formation of planets in the HZ. Ice lines are also in relatively inner 
regions at ~ 0.3AU. Due to the small orbital radii, icy protoplanets accrete rapidly and 
undergo type-I migration before disk depletion. The rapid orbital decay, the proximity 
of the disk inner edge, and large amount of inflow of icy protoplanets are characteristic 
in planetary accretion in terrestrial planet regions around M dwarfs. In the case of full 
efficiency of type-I migration predicted by the linear theory, we found that protoplanets 
that migrate to the vicinity of the host star undergo close scatterings and collisions, and 
4 to 6 planets eventually remain in mutual mean motion resonances and their orbits 
have small eccentricities (< 0.01) and they are stable both before and after disk gas 
decays. In the case of slow migration, the resonant capture is so efficient that densely- 
packed ~ 40 small protoplanets remain in mutual mean motion resonances. In this case, 
they start orbit crossing, after the disk gas decays and eccentricity damping due to tidal 
interaction with gas is no more effective. Through merging of the protoplanets, several 
planets in widely-separated non-resonant orbits with relatively large eccentricities (~ 
0.05) are formed. Thus, the final orbital configurations (separations, resonant or non- 
resonant, eccentricity, distribution) of the terrestrial planets around M dwarfs sensitively 
depend on strength of type-I migration. We also found that large amount of water-ice 
is delivered by type-I migration from outer regions and final planets near the inner disk 
edge around M dwarfs are generally abundant in water-ice except for the innermost 
one that is shielded by the outer planets, unless type-I migration speed is reduced by a 
factor of more than 100 from that predicted by the linear theory. 
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1. INTRODUCTION 



Over 300 extrasolar planets have been discovered. Target stars for exoplanet search were 
mostly solar-type stars (F, G, K dwarfs), although M dwarfs make up 70-80% of all stars in the 
galactic disk. The low luminosity of M dwarfs is disadvantageous for high-dispersion spectroscopic 
observation, so that radial velocity surveys have not discovered large number of planets around 
M dwarfs. However, as improvement of spectroscopic observations, ground-based radial velocity 
surveys are revealing planetary systems around M dwarfs. Due to the low luminosity of M dwarfs, 
habitable zones (HZ), in which a planet with su fficient amount of atmosphere can sustain liquid 
water on its surface, are close to the host stars (iKasting et al. |[l993). The proximity of the HZs 
to the host stars allow for detection of planets in HZs by current radial velocity observation. In 
fact, two planets with minimum masses below lOM^ were discovered near the HZ in a triple planet 
system around an M star, Gliese 581, with stellar mass = O.SIM© (jUdrv et al.l 120071 ). The 
habitability of these planets (Gl 581c, d) is vigorously under discussion theoretically. In addition, 
gravitational microlensing survey is suited for detection of M dwarf planets, because its detection 
efficiency is independent of stellar luminosity. Most of the planets detected by microlensing are 
orbiting M dwarfs. Recent radia l velocity and ni i crolensing observati ons show that Jupiter-mass 



are rather abundant (e.g.. iBeaulieu et al.ll2006l ). compared with solar-type stars. 



gas giants are generally ra re fe.g.JEndl etal.ll200d . I Johnson et al.ll2007l ). but Neptune- mass planets 



GJ 436b is the planet that was discovered first among Neptune-mass planets (jButler et al 



20041 ). Transit observations revealed a planet's radius, and its combination with radial velocity 
measurements permits a determination of the planet' s density. The evaluated internal den sity 
suggests that GJ 436b can be composed mainly of ice (jGillon et al.l 120071 . iDeming et al.l 120071 ). in 
spite of proximity to the host star. On-going and upcoming transit surveys using space telescopes 
such as Corot, Kepler and TESS, besides ground-based transit surveys, are expected to reveal lower 
mass exoplanets around M dwarfs. 



"Core accretion" model (e.g., iHavashi et al.lll985l ) naturally accounts for the low abundance 



of gas giants around M dwarfs, because observationally inferred low disk mass around M dwarfs 
inhib its formation of cores large enough for runaway gas accretion (ILaughlin et al.ll2004bl : llda &: Lin 
20051 ). The relatively high abundance of Neptune-mass planets could be acc ounted for by tru ncation 
of gas accretion at smaller plane tary mass due to lower disk temperature (jlda &: Linll2005l ). Thus, 
the Monte Carlo calculation by llda &: LinI (120051 ) (at least qualitatively) explained the observed 
properties in the population of gas giants and Neptune-mass planets around M dwarfs. However, 
they did not predict abundance of habitable planets, because it is affected by type-I migration 
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and detailed orbital configurat ions of close-in planets, which were not taken into account in their 
calculation. Ilda &: LinI (|2008al ) included the effect of type-I migration in the similar calculation, 
but treatment of close-in planets was still too simple to discuss the abundance of habitable planets 
at ~ 0.1 AU around M dwarfs. 

A'^-body simulation is an efficient tool to address this issue. Since physical sizes of planetesimals 
occupy a larger fraction of their Hill radii in the terrestrial planet regions (~ 0.1 AU) around M 
dwarfs than around solar-type stars, strong gravitational scattering is suppressed, which reduces 
computational cost. Moreover, we can neglect perturbations from gas giants because they are rare 
in the M-dwarf planetary systems, which also makes iV-body simulation simple. 



Raymond et al.l (|2007l ) performed a first iV-body simulation of terrestrial planet formation from 



planetary embryos around low-mass stars. They found that the planets in a HZ may be too small to 
retain ocean because they assumed that disk surface density is proportional to the stellar m ass and 
the disk model for IMq is the minimum mass solar nebula (MMSN) model (|Havashilll98ll ). Under 



this assumption, the isolation mass of the planets is proportional to stellar mass (section 2.3). In 
HZs, icy grains do not condense in the protoplanetary disk in which gas pressure is much smaller 
than the planetary atmospheric pressure. One of available sourc es for the water on the planets is 
delivery of icy planetesimals from the regions beyond the ice line (iMorbidelli et al.ll200Cll ). although 
the possibility of forming H2O through chemical i nteraction between the planetary magma ocean 
and primitive H2 atmospher e is also pointed out (Ikoma &: Gendal l2006l ) . Assuming the delivery 
hypothesis of origin of H2O, [Raymond et al.l (j2007l ) suggested that the planets in HZs around M 
dwarf stars are likely to be dry, sinc e radial mixing and therefore water delivery are inefficient in 
the lower-mass disks. iLissaueri (j2007l ) also pointed out the possible lack of large volatile inventories 
due to the large collision speeds of impacting comets and substantial mass loss of volatiles due to 
high activity and luminosities of young M dwarfs. 



Although the simulations by [Raymond et al.l (120071 ) provided important insights into proba- 
bility of habitable planets around M dwarfs, they neglected the effects of protoplanetary disk gas. 
Since in such inner regions, gas density is so high that migration due to gas drag and tidal inter- 
action with a gas disk are efficient (see section 2.3) and accretion timescale of terrestrial planets 
would be much shorter than disk lifetime (see section 2.3), the effects of disk gas play important 
roles in the accretion of planets in HZs around M dwarfs and water delivery to them. We will point 
out in section 3.5 that the planets in HZs are rather composed mainly of water-ice unless type-I 
migration speed is reduced by a factor of more than 100 fr om the linear theory or the migrating 



protoplanets are tr apped at an inner boundary of dead zone ([Kretke &: Lin 
Kretke et allboOfll ). 



2007 



Ida Linll2008bl . 



With type-I migration, the proximity of the HZ to the inner disk edge would play an important 
role for final configuration of pla nets in HZs, because type-I rn igration stops at the disk edge and 
planets would accumulate there. iTerquem &: Papaloizoul (|2007l ) performed iV-body simulations of 
protoplanets undergoing type-I migration around solar-type stars. Their important finding is that 
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the migrating protoplanets originally formed at ~ lAU interact with the preceding protoplanets 
near the disk inner edge and finally two to five close-in planets remain in mutual mean-motion 
resonances. They found that the resonant configuration is maintained even after disk gas is removed 
or tidal interaction with the star is added. They only carr ied out simulations with type-I migration 
with full efficiency that is predicted by the linear theory (jWardlll986l . iTanaka et al.ll2002l ). We will 
show that the final orbital configuration is sensitively dependent of migration speed. 

We thereby carry out A'^-body simulations including the effects of damping of orbital eccentric 



ity, iri clina tion and semimajor axis (type -I migration) due to disk gas. Although iRavmond et al 



(|2007l ) and lTerquem &: Papaloizoul (|2007l ) start their simulations from planetary embryos that have 
already grown to their isolation masses, our simulation starts from many small planetesimals, tak- 
ing fully into account their gravitational interactions. Our calculations also cover a broad range of 
orbital radius from the disk inner edge to beyond the ice line. Since there is still uncertainty in the 
type-I migration speed, we perform both simulations with and without type-I migration. 

In section [2l we describe the disk model, the formulas of forces for aerodynamic and gravita- 
tional gas drag, and calculation methods. The results of iV-body simulations of planetary accretion 
are shown in section [3l In section 13.41 we discuss water delivery to inner planets around M dwarfs. 
Section H] is devoted to the conclusion and discussion. 



2. MODEL AND CALCULATION METHODS 

Here, we consider an M5V type star with mass M^, = {).2Mq and luminosity L^, ~ O.OILq, 
using a mass-luminosity relation (L* oc M^). This relation roughly fits observational data for 



a ran ge of 0.1 to 1 Mq stars in main-sequence stages (e.g., iHabets &: Heintzd Il98ll . IScalo et al 



20071 ) . Although pre-main-sequence stages are relatively long for the low mass stars and luminosity 
is relatively high during the pre-main-sequence stages, we will use HZs determined by the main- 
sequence radiation, because our main purpose is to clarify the dynamics and accretion process 
among protoplanets that have migrated to the regions near disk inner edges. The simulation with 
evolving HZs and ice lines due to the luminosity evolution is left to a future work. 

We integrate the orbits of planetesimals, taking into account their merging by direct collisions, 
their gravitational interactions, aerodynamic gas drag, "gravitational drag" (damping of orbital 
eccentricity and inclination due to tidal interaction with a gas disk), and type-I migration. The 
models for surface densities of disk gas and an initial planetesimal swarm are explained in section 
2.1. Basic equations for orbital integration and initial conditions are presented in section 2.2. 
Although the detailed expressions for the aerodynamic gas and gravitational drag forces are given 
in Appendixes A and B, we summarize their characteristic timescales in section 2.3, as well as 
the timescale of planetesimal accretion, which are useful to understand the results of the A^-body 
simulations. 
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2.1. Disk Model 



Following llda &; Liij (j2004l ) , we scale the gas surface density of disks as 

-3/2 



2400/, 



lAU 



g cm 



(1) 



where fg is a scaling factor; is 1.4 times of those in the MMSN model if fg = 1. Because current 
observations cannot strictly constrain the radial gradient of Tig, we here assume fg is constant with 
r except for the inner edge. For disks around stars with ~ ^Mq, the observationally inferred 
averaged val ue of fg is ^ 1 a lthough the values have dispersion of two orders of magnitude (see 
discussion in llda Linll2004l ). We here adopt fg = l for the disks around the M* = O.2M0 star. 



Although averaged fg may be several times smaller for these stars, = 1 is within the two orders 
of magnitude dispersion. The relatively large value of fg is to study possible formation of habitable 
planets, which are large enough to retain water on their surface, around low mass stars. In section 
3.5, we will discuss how the results are affected if less massive disks {fg ~ 0.2), which may be 
averaged disks around these stars, are considered. 



We assume the temperature distribution of an optically thin disk (|Havashilll981 



T = 2.8 X 10^/3 
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Corresponding sound velocity Cg and disk scale height h are 
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In the current paradigm, T Tauri disks are truncated by the stellar magnetosphere within the 
corotation radius, with materials accreting along magnetic field lines onto high-latitude regions of 
the star. The corotation radius rcorot is the radius at which the Keplerian orbital period in the disk 
equals the stellar rotation period (P), 



^'corot = 0.04 



3days 



(5) 



0' 



Here, we set the location of disk inner edge at 0.05 AU. From Eq. (UD, the scale height of a disk at 
0.05 AU around a star with M,, = 0.2Mq is ~ 1.4 x 10""^ AU. We assume that the gas surface density 
of the disk (equivalently fg) smoothly vanishes at 0.05 AU with a hyperbolic tangent function with 
the width of 10~^ AU, which is comparable to h. When planets enter the inner cavity, they do not 
feel gas drag any more and their inward migration ceases because the drag and migration rates are 
proportional to disk gas surface density. In some runs, w e reverse the direct i on of the migratio n 
near the edge according to positive pressure gradient there (jTanaka et al.ll2002l : lMasset et al.ll2006l ). 
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According to gas surface density given by Eq. ([T]) , we scale the surface density of a planetesimal 
disk with a scahng factor, f^, as 



Srf = lOr/ice/d 



lAU 



-3/2 



gem 



(6) 



For solar metallicity, fa = fg, so we use = 1. The factor r/ice expresses the increase of solid 
materials due to ice condensation outside the "ice line." The location of the ice line is determined 
in such a way that the disk temperature equals the ice condensation temperature. Assuming the 
condensation temperature of 170 K, the location of the ice line is derived from Eq. ([2]): 



2.7 



Since we consider stars with L 

at r > rice 



O.OIM0, we set rice 



AU. 



(7) 



0.3AU. In the MMSN model, r/ice = 4.2 



Pollack et al.l derived 

^icc — 3. In the standard set of our simulations, we adopt 



as 



[r < 0.3AU] 
[0.3AU < rl. 



(8) 



But, in some runs, we adopt a higher value for rice (see section 3.3). 



2.2. Orbital Integration and Initial Conditions 



We integrate the orbits of planetesimals with 4th-orde r Hermite scheme (jMakino &: Aarseth 



1992) and hierarchical individual timestep (iMakind Il99ll ) . The basic equations of motions of 



particle k at in heliocentric coordinates are 



~l~-^aero ~l" -f^damp ~l~ -f^migj (9) 

where k,j = 1,2,..., the first term is gravitational force of the central star and the second term is 
mutual gravity between the bodies. We calculate the self gravity directly summing up interactions 
of all pairs on the special-purpose computer for TV-body simulation, GRAPE-6. -Faeroj -P'damp and 
Fmig are specific forces due to aerodynamic gas drag, gravitational drag, and type-I migration, the 
detailed expressions of which are described in Appendixes [X] and [Bj We neglect the indirect term 
since the total mass of the planetesimals is ~ 10~^ times the mass of the central star. 

When physical radii of two bodies overlap, perfect accretion is assumed. After the collision, 
a new body is created, conserving total mass and momentum of the two colliding bodies. The 
physical radius of a body is determined by its mass M and internal density pp as 

/ 3 M\i/3 , ^ 
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We adopt a realistic value 3 gem ^ for pp. 

Initially 5,000 planetesimals are placed between 0.05 AU (disk inner edge) and 0.4 AU. To 
study accretion process of terrestrial planets inside the ice line in more detail, we initially set 
more bodies inside the ice line (3898 bodies with mass 2.3 x 10^4 g) than outside it (1102 bodies 
with mass 6.5 x 10^^ g) in the nominal simulations. The initial velocity dispersion of the bodies 
is set to be their escape velocity. The corresponding in itial eccentricity and inclination are given 



by Vesc = Ve^ + i'^VK with e = 2i ( Ida &: Makindll992l ). Table [T] shows simulation conditions for 



individual runs and final results. In the first set (set A) of runs (runAl - runA4), we include all 
the damping of e, i and a due to aerodynamic and gravitational drag and type-I migration, while 
type-I migration is neglected in the second set (setB: runBl - runB4). Because coagulation to 
final planets may be a stochastic process, for each set we performed 4 runs with different seeds for 
random numbers to generate initial angular distribution of planetesimals. 

In addition to these sets, we also carried out runs with reversed torque of type-I migration 
near the disk inner edge (setC: runCl - runC2) and runs with rjice = 14 at r > rice (setD: runDl - 
runD4 with type-I migration, setE: runEl - runE4 without type-I migration), which are described 
in section 3.3. 



2.3. Characteristic Timescales 

To understand the results of TV-body simulations, we here summarize the timescales of gas 
and gravitational drag, type-I migration and planetesimal accretion. The characteristic damping 
timescale for orbital eccentricity (e) and inclination (i) of an isolated p lanetesimal due to aerody- 



namic gas drag is (for more detailed expressions, see lAdachi et al.lll976l ) 



laero ~ „ _ O.^ X iU ^ q j ^Mfs J \0.2Mq J 



Pp 



N 2/3 / N 13/4 

^) [ih) years, (11) 

where Faero is specific drag force acting on the planetesimal (eq. |Alj ) and a, M and pp are its 
semimajor axis, mass and bulk density, respectively. The relative velocity between the planetesimals 
and disk gas, Au, is given by ~ ((5/8)e2 + (1/2)^2 + ri'^)'^/^VK, where vk is Keplerian velocity. The 
velocity of disk gas Ugas is smaller than Kepler velocity vk by a fraction (eq. |A2j ) . 

JQ. 

The timescale for damping of semimajor axis (a) is 



,^^^2,8xl0-'(— j . (12) 
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At a ~ O.IAU, taero,a ~ 1-5 X 10^ years even for a Mars-mass planet (M ~ O.IM®), so that 
gas drag cannot be neglected even for protoplanets. For small planetesimals, the gas drag is more 
important. 

A planet gravitationally pertur bs the disk gas and excites de nsity waves. The waves damp 



the eccentric ity, the inclination (e.g. , Ward 19931 . Artvmowicz 1993 ) and the semimajor axis of the 



planet (e.g.. IColdreich fc Tremain] Il98d . IwardI Il98d l. The detailed expressions of gravitational 
gas drag forces -Fdamp and -Fmig are given in Appendix [Bj Orbital eccentricities and inclinations 
are damped by both torques from inner and outer disks in a similar way to dynamical friction from 
planetesimals. The damping timescales are (jTanaka h WardI 12004 ) 
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On the other hand, secular inward migration due to tidal interaction with disk gas, that is known 
as "type-I migration" is caused by torque imbalance. The migration timescale is 
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(16) 
(17) 



assuming the gas surface density is proportional to a Here we used the MMSN density profile, 



Tig (X a , to derive Eq. (I17p . In runCl and runC2, the effect of negative q at the inner edge, 
which reverses the tidal torque, is taken into account. 

Since type-I migration is resulted in by torque imbalance, non-linear effects could change 
migration rate significantly. Therefore we also perform simulations in which bodies are free of 
type-I migration (runBl - runB4). Note that even in this set, the migration induced by damping 
of e due to the gravitational drag exists (section 3.1.2) and the runs in this set are comparable to 
the runs with type-I migration of ~ 100 times reduced speed. 

Growth rate of a protoplanet with mass M and physical radius rp due to accretion of small 
planetesimals is estimated simply by 
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K, 



(18) 



where Uran — {e^ + i'^f/'^VK- With more detailed formula and t^ran that is determined by balance 
between scattering of the planetesimals b y the protoplanet and gas drag to them, the accretion 
timescale is expressed as (jida h Linll2004l ) 
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The growth timescale of terrestrial planets around M dwarfs is shorter than that around solar- 
type stars because of small a. Hence, it is reasonable to assume that full amount of disk gas exists 
during the terrestrial planet formation because depletion timescale of disks around M dwarfs may 
be ~ 10^ — lO''' years or more. On the other hand, the damping timescales are short in inner regions 
where the gas density is high. As stated in section 1, the effects of disk gas cannot be neglected 
when we consider accretion of planets in HZs around M dwarfs. 

The isolation mass, which is the mass of all the solid materials in a feeding zone of the proto- 
planets, is 

nOQ 3/2,3/2/ a \3/4/ Aa \3/2/ \-i/2 

where Aa is width of the feeding zone. The mutual Hill radius th is defined by 

_ / Ml + M2 \V3 Miai +M2a2 

where Mi and M2 are masses of interacting bodies that we are concerned with and ai and 02 are 
their semimajor axes. In Eq. (I20p . Mi = M2 = Miso is assumed. The value of Aa = 7.5rH is a 
typical value obtaine d by A^-body simulation at a ~ 0.1 AU, which is slightly smaller than the value 



obtained at ~ 1 AU (jKokubo &: Idalll998l . |2002| ) . The isolation mass is the maximum mass that the 
planet can acquire through "runaway/oligarchic" growth (accretion of planetesimals) before onset 
of orbit crossing and coagulation among the isolated protoplanets. 

The migration timescale for the protoplanet with mass Mjso is 

3.x:0*Cr//''/-(^)'"(^)"%ea.. (22) 



Its accretion timescale is 



Wo - U.70 X 10 r/i^, ^^^^j U.2MqJ Ugcm-3j y"^''' ^^^^ 

At a < 0.3AU that we are concerned with in this paper, both tmig,iso and tacc,iso are significantly 
shorter than disk lifetime of ~ 10^ — 10'' years for /g ~ /d ~ 1. This is also the case even for 
small- mass disks with t he su rface density 5 times smaller than that of the MMSN {fd = fg — 0.15), 



which [Raymond et al.l (120071 ) considered. Thus, type-I migration cannot be neglected even for the 
small-mass disks, unless non-linear effects slow down or halt the migration significantly (section 
3.5). 
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3. RESULTS 

3.1. Runaway /oligarchic Growth before Disk Gas Depletion 

3.1.1. Case with Type- 1 Migration 

The result for a typical run including the effect of type-I migration (runAl) is shown in Figured) 
In the snapshots, the sizes of circles are proportional to the physical radii of the bodies, and Tk is 
the Kepler time at 0.1 AU around a 0.2Mq star, which is ~ 0.071 year. Since the accretion timescale 
depends strongly on semimajor axis a (Eq. [19]) and it is shorter for smaller a, accumulation of 
planetesimals proceeds in an inside-out manner. Equation (jl9p for M ~ M^^o ~ O.O4M0 at 0.1 AU 
is 10^ yrs = 1.4 x 10^ Tk, which agrees with the result in Fig. [TJ 

According to growth of a protoplanet, the type-I migration timescale becomes shorter, while 
the accretion timescale becomes longer. Thereby, the protoplanet starts migration when its mass 
exceeds the critical mass beyond which tmig < tacc- From Eqs. (fT7)) and (fT9|) . the critical mass is 

crit.mig /ice ./d Jg VlAU/ VO.2M0/ VO.OILq/ ® ^ ' 

Because at ~ O.IAU this value is close to the isolation mass, protoplanets start migration after 
they accrete most of planetesimals in their feeding zone. On the other hand, in outer regions, 
since Mcrit,mig < -Miso, protoplanets start migration leaving behind large amount of planetesimals. 
Most of the initial mass is finally concentrated near the disk inner edge after 10^ Tk, leaving few 
planetesimals in outer regions. Although eccentricities of the bodies are excited by close encounters 
with neighbor bodies during the growth stage, final eccentricities of the planets are kept small 
(< 0.01) due to gravitational drag. Resonant effect can raise the eccentricities, however, gas drag 
damps them significantly because of the high gas surface density. 

Orbital evolution of runAl is shown in Fig. [2j In this figure, the most massive 30 planets at 
each time are plotted as circles (there are many other small bodies). Accretion takes place more 
rapidly in inner region than in outer region at the beginning of the simulation. The protoplanets 
with M > Merit, mig Undergo inward type-I migration. Since tmig is shorter at smaller orbital 
radius, the migration is accelerated until the protoplanets reach the disk inner edge. The firstly 
migrated protoplanets interact with each other. After some merging, the survived protoplanets are 
captured in mutual mean motion resonances. After that, successively migrated planets interact 
with outermost planets that have accumulated near the disk edge. The interaction mostly results 
in merging of the planets and the merger is again captured in a mean motion resonance. Note that 
the number of remaining planets near the edge is almost constant (4-6 planets) during this phase, 
although protoplanets that formed in outer regions migrate to interact with the close-in planets 
one after another. In this run, after ~ 5 x 10^ Tk, any more protoplanets which are large enough 
to affect the inner planets do not approach to the inner planetary system, so that this run ends up 
with a stable configuration consisting of six planets with M > O.OIM® near the disk inner edge. 
Most of the final planets are pushed into the cavity at a < 0.05AU, in which type-I migration is no 
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more effective, by the outer planets that keep loosing angular momentum by type-I migration. The 
largest planet in the final state is the third innermost planet, the mass of which is O.63M0. Note 
that this value is more than 10 times of the values of Mcrit,mig (~ 0.01 — O.OGM® at 0.05-0.4AU). 
This implies that coagulation near the disk edge is so efficient. The value of O.GSMq is also more 
than 25 times larger than Miso (~ 0.024:M^) at 0.05AU. Without migration, such large planets 
cannot accrete at ~ 0.05AU. 

All the final six planets are trapped in first-order commensurability (mean-motion resonances) 
with the planets which lie next to the planets. For example, the innermost pair (a pair of the 
innermost and the second innermost planets) has 6:5 commensurability, and the second pair (a pair 
of the second and the third innermost planet s) has 5:4 commensurability. The perturbations during 



one passage rapidly increase for Aa < 5rH (jldalll99d ). Since resonant trapping at Aa > 5rH may 



not be resistant against perturbations from other bodies other than the pair, it is expected that 
the pair may be eventually trapped at mean-motion resonances close to Aa ~ 5rH. The obtained 
orbital separations of the i-th pairs {i = 1,2,. .,5) are Aa ~ 7.3rH, S.Srn, 6.8rH, 5.4rH and 8.7rH, 
respectively. As will be shown in section 3.2, these final orbital configuration is stable even if disk 
gas is removed. 

All the runs in set A with the effect of type-I migration end in very similar results: Most of the 
protoplanets (large planetesimals) pile up near the inner disk edge, and in final state, most of the 
planets are captured in mean motion resonances after coagulation and close scattering among the 
protoplanets. The average numbe r of the final planets in r unAl - runA4 is A'^ = 5.0 it 0.71, which 



is comparable to that obtained by iTerquem &: Papaloizoul (120071 ) . The average mass of the largest 



planet is M^ax = 0.82 ± 0.17M©, which is significantly larger than Miso &t 0.05AU and -Merit, mig 
at 0.05-0.4AU. 



3.1.2. Case without Type-I Migration 

A typical result excluding the effect of type-I migration in setB is shown in Figure [3l The figure 
shows snapshots of the system on the a — e planes for runBl. Although the effect of type-I migration 
is not included, planets tend to migrate inward. The inward migration is induced by eccentricity 
damping by tidal interaction with disk gas. Since the angular momentum (L = y^GM^,a(l — e^)) 
is almost conserved, damping of eccentricity yields damping of semimajor axis. The damping 
timescale of semimajor axis is roughly given by 



0.7xl06f-M777T7 (t^t) TTTT 77777^ years, (26) 



» VO.Ol/ VMe/ VIAU/ VO.2M0/ Vo.OILq 

where Eq. was used. The values of e are typically ~ 0.01 in our simulations, but they change 
with time and depend on locations. Compared to Eq. (jl3|) . we find that the semimajor axis 
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damping induced from the eccentricity damping is more effective than aerodynamic gas drag for 
M ~ 0.01 — IM® , even if the uncertainty in the values of e is taken into account. Since this timescale 
is about 100 times longer than tmig (Eq. [T7]) at a ~ 0.1 AU, this calculation is equivalent to the 
case with 100 times reduced type-I migration speed. Nevertheless, this migration timescale is of 
the order of lO^yrs for a planet of O.IM^ at 0.1 AU around a star of O.OILq, which is considerably 
shorter than disk lifetime. This effect plays an important role in planetary formation in HZs 
around M dwarfs, although it can be neglected in HZs (a ^ lAU) around G dwarfs. This migration 



phenomenon can be seen in previous works including tidal e-damping with gas disk (jOgihara et al 



20071 ). In Fig. [31 the planets have a tendency to move inward. However, several planets remain in 
outside regions whereas in runAl few planets remained outside region. The final eccentricities are 
damped down to 0.01 as in the previous set. 

Figure [4] is the same plots as Fig. [2] for runBl. In this run, the resonant capture is so effective 
because of the slow migration that about 45 planets with mass larger than 0.01 are lined up in 
resonances with little coagulation among them. On the other hand, in runAl, migrated planets were 
not readily trapped at a first encounter and settled to resonant configurations after close scattering 
and coagulation. In runBl, the largest planet is located at 0.21 AU, the mass of which is 0.21 M^. 
Almost all the planets are captured in mean motion resonances, however, they tend to have closer 
commensurabilities, such as 14:13, than those for runAl. In the case with type-I migration (setA), 
successively migrated planets interact with outermost planets, leading to merger of planets. In this 
case, however, perturbations from many other migrated planets tend to break resonant capture at 
distant separation, resulting in capture in closer (stronger) mean motion resonances. 

All the simulations in this set (runBl - runB4) exhibit qualitatively the same results. The 
average number of final planets is as large as = 40 it 3.3, which is about one order larger than 
that in setA. The average mass of the largest planet is Mmax = 0.20 it O.O33M0. 



The critical planet mass beyond which tdamp.a < ^acc is 

-3/2^ a \-V4/ \3/8 

f©/ Vo.OILq- 



3/4.3/4,_9/20/ e X-3/2/ a X -21/40 /M^x -1/4 .^.X 3/8 



This value is greater than the isolation mass (Eq. [20]) at a < 2AU, so that all the protoplanets 
in the simulations in this set grow up to the isolation mass before onset of migration. Hence, the 
maximum mass Mmax is roughly equal to Afiso, which is consistent with the value of Mmax obtained 
in our results. 

The orbital separations are ~ 5 — Grn, which means that the planets are more packed than 
in setA. The larger number and smaller separations of the final planets than those in setA suggest 
that the final configuration could be unstable. As will be shown below, in setB, the final planets 
start orbit crossing immediately after the removal of disk gas. 
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3.2. Stability after Disk Gas Depletion 



All the runs in setA with type-I migration and setB without it end with multiple planets. 
Because eccentricity damping due to gravitational drag i s so strong that the sy stems of these 



multiple planets are orbitally stable (jlwasaki et al. 



2002 



Kominami &: Idal |2002| ) , although the 



minimum orbital separations in the systems are rather small ('- 
dissipate on timescales less than 10'' years. 



5 — Trn). However, disk gas should 



In the gas-free case, it is predicted that these p lanets may become unstable on timescales 
of tcross ~ IO^Tk ~ 10^ years ( Chambers et al.lll996l ). Note, however, that these timescales are 



for non-resonant planets. In our results, the final multiple planets are usually captured in mean- 
motion resonances that generally stabilize the systems. Actually, iTerquem &: Papaloizoul (|2007l ) 
found that the final multiple planets, which correspond to our results in runAl - runA4, are stable 
on timescales much longer than icross even after removal of disk gas. We show similar results for 
runAl - runA4 below, but also show that the removal of disk gas makes the systems unstable for 
runBl - runB4. 

To examine long term stability, we take out the planets with masses larger than O.OIM^ from 
the final state of the runs, and integrate their orbits, neglecting many other smaller-mass planets 
for saving computational cost. As expected, we found that the system of runAl is stable until 
2 X IO^Tk under the disk gas damping. To study the stability after the disk gas removal, we re- 
start the calculation from the final state of runAl at 5 x IO^Tk without the disk gas damping. We 
adiabatically adjust the orbital configuration to the gas- free condition, by decreasing fg as 

t - 5 X IO^Tk ' 



fg = exp 



lO^years 



(28) 



The decay timescale of 10^years(~ 1.4 x IO^Tk) is long enough for the adjustment. 



At the end of runAl (Fig. six planets are remain with separations of 5 — Orn- The orbital 
evolution after the disk gas removal of runAl is shown in Fig. [5] (the left panel). The figure shows 
that the eccentricities of the planets are kept less than 0.01 and the planets remain stable even after 
the gas removal. Although the minimum orbital separation is about 5rH, the resonant configuration 
stabilizes the system. 

The other runs (runA2 - runA4) in this set with type-I migration also show similar results: 
Even after the disk gas removal, orbital separations hardly change (the average orbital separation 
is Aa = 9.5 it 0.97rH) and most of the planets keep their commensurate relationships until the end 
of simula tions. The average orbital ecc entricity is e = 0.0086 ± 0.0061. This is consistent with the 
result by iTerquem &: Papaloizoul (120071 ). 



The stability for setB without type-I migration is completely different. At the end of runBl 
(Fig. S]), 45 planets with masses larger than O.OIM® remain with the orbital separations of the 
planets are 5 — 6rH. Figure [5] (the right panel) shows the semimajor axis evolution of runBl after 
the disk gas removal. Soon after the gas removal, the eccentricities are pumped up and the planets 
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start orbit crossing. Note, however, that the planets do not exhibit global orb ital crossing that 
the p lanets at ~ lAU around solar-type stars exhibit after disk gas removal (e.g.. iKominami &: Ida 



2002 ). In terrestrial planet regions around M dwarfs (~ O.IAU), physical radii of the planets 
Tp relative to their Hill radii rn are larger than that in ~ lAU (rp/rn oc M^^^a"^). Then, the 
eccentricities are pumped less highly and furthermore merging proceeds before the eccentricities are 
fully excited. Thus, the planets collide with only neighboring planets. Finally, nine planets with 
moderate eccentricities 0.08) are formed. The mass of the largest planet is 0.65M^. All the 
commensurabilities are lost in the course of close encounters and the final planets are not trapped 
in mean motion resonances at all. The orbital separations are ~ 20rH, which are large enough 
to be dynamically isolated from each other in the non-resonant configurations. Since relatively 
many planets (~ 40) are formed with small orbital separations, all the runs in setB without type-I 
migration exhibit orbital crossing and merging of the planets after the disk gas removal, resulting in 
~ 10 — 20 planets with the average eccentricity of 0.055 it 0.020 and the average orbital separation 
of 19 lb 2.2rH that lost commensurate relationships. The average mass of the largest planet is 
-^max = 0.50 lb 0.097-/Vf^. The non-resonant orbital configurations with wider separations and 
larger eccentricities are the characteristic to setB without type-I migration (with 100 times reduced 
migration efficiency). The semimajor axes of the planets are not concentrated to the regions near the 
disk inner edge, in contrast to the results with type-I migration. Thus, type-I migration efficiency 
in inner disk regions regulates orbital configurations of close-in terrestrial planets. 



3.3. Dependence on Boundary and Initial Conditions 



Since we are concerned with planetesimal accretion near the disk inner edge, we study the 
effects of general relativity that is effective in the proximity of the host star and reserved type-I 
migration torque that occurs near the edge. We re-started the stability calculation in gas-free 
condition for the runAl - runA4, incorporating the relativistic effect directly into orbital integra- 
tion. The detailed expression of post-Newtonian gravitational force from the host star is given in 
Appendix O Although the relativistic effect causes the precession of the perihelion of short-period 
planets, we found that the resona nt relationships are not change d by the relativistic effect and the 
systems stayed in a stable state. iTer quern fc Papaloizoul (j2007l ) found that the tidal dissipation 
does not affect the stability as well as the relativity. 

It is argued that inward protoplanet migration can be halted before reaching the inner cavity, 
because the tidal to rque from the disk is reversed due to inverse pressure gradient near the disk edge 
(jMasset et al.ll2006l ). Because the planets in the inverse torque regions gain angular momentum 
from the disk gas, the inwardly migrating planets in outer regions cannot push the inner planets into 
the cavity, before the depletion of disk gas. Substituting the component of the gas surface density 
gradient at the inner edge in our model into —q in Eg. (I16D. we performed simulations in runCl and 
runC2 (This effect was also investigated in iTer quern fc Papaloizoul (j2007l )). The orbital evolution 
for runCl is shown in Fig. [6l Although qualitative evolution is almost the same as the result 
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without the reversed torque, the number of final planets is fewer than that obtained in runAl 
- runA4. Because the innermost planet cannot penetrate into the cavity, the orbital separation 
between the innermost planet and the second innermost one is small (~ 3.8rH). However, the right 
panel of Fig. [6] shows that the planets remain orbitally stable after the disk gas removal, because 
they keep being trapped in 5:4 and 9:8 resonances. 

Around M dwarfs, the ice line is so close that significant amount of icy protoplanets are quickly 
formed and migrate to inner regions. The migrated icy protoplanets affect accretion of planets in 
terrestrial planet regions, in particular, they regulate the mass of the largest planet in final state. 
So far, we have adopted the ice condens ation factor as r]\c.p, = 3 in Eg . ([6]). However, more enhanced 



r/ice was proposed by several authors. IStevenson &: Lunind (jl988l ) proposed that sublimation of 



icy grains that have migrated to inside the ice line and the diffusion of the water vapor enhances 
the surface density of icy materials near the ice line up to jyice ~ 75. Recent more detailed study 



(ICiesla Sz Cuzzil l2006l ) showed ??ice ~ 10. To highlight the effect of migrated icy protoplanets, we 



performed calculations with 



r 1 [r < 0.3AU] 

I 14 [0.3AU < rl. ^ ^ 



"^''^ ' 14 [0.3AU<r] 

The result including the effect of type-I migration (runDl) is shown in Figure [71 It is qual- 
itatively similar to runAl with r/ice = 3 at r > 0.3AU: several planets are formed near the disk 
inner edge trapped in mutual mean motion resonances before the disk removal. Even if disk gas 
is removed, the final planets are stable. We also did three additional runs with the same setting 
(runD2 - runD4) and they show similar results. However, the average mass of the largest planet is 
-^max = 3.6 lb 0.29Mq), which is significantly larger than that with r/ice = 3 (setA). This difference 
means that the planets mostly consist of migrated icy protoplanets. Because orbital separations 
must be greater than ~ Srjj for the planets to be stable and rjj is proportional to M^/^, the average 
number of final planets {N = 4.0 it 0.71) is smaller than that in setA. 

We also performed runs with rjice = 14 that did not include type-I migration (setE; runEl - 
runE4). The result for runEl is shown in Fig. [HI Before the disk gas removal, the results in setE are 
similar to those in setB with r/ice = 3, except for the final maximum mass (Mmax = O.91ibO.O64M0) 
and number (N = 27 it 2.3). In setB, resonant trapping was so efficient that Mmax was no other 
than Merit, mig- In setE, the inner protoplanets cannot halt the migrated protoplanets from outer 
regions, total mass of which is one order larger than the total mass of inner planets, and mergers 
and rearrangement occur in the inner regions. This results in the larger Mmax and smaller N in 
setE than in setB. 

So far, we have been using Eg oc r~^'^. We also did several simulations with less steep 
radial gradient, oc r~^'^. Because of weaker dependence of type-I migration timescale on r 
corresponding to the weaker r-depende nce of S^,, we found t hat a few protoplanets that are trapped 



by each other migrate together, which [McNeil et al.l (|2005[ ) called a "convoy." But, we also found 



that this feature does not affect the final orbital configurations of close-in planets before disk gas 



removal and their stability after the disk gas removal. The final orbital configurations depend 
on only type-I migration speed around O.IAU, because it determines the efficiency of resonant 
trapping. 



3.4. Composition and Habitability 



As stated in section 1, delivery of icy planetesimals fro m the regions beyond the ice line is 
one of likely sources for the H2Q-water o n plan etary surface (jMorbidelh et ahlboool . iRobertlboOlh . 
Assuming this scenario, iRavmond et al.l (|2007l ) suggested through iV-body simulation neglecting 
type-I migration that the planets in HZs around M dwarf stars are likely to be dry, since radial 
mixing is inefficient in the lower-mass disks. 

Figure [9] shows H20-water mass fraction of the final planets in runAl (the left panel) and 
runBl (the right panel), using the following simple prescription for components of planetesimals 
that originated at r: 



water 



M 



I 0. 



[r < Tice] 

67 [r > rice]. 



(30) 



Note that the initial mass beyond the ice line make up 41% of total mass in our calculation range, 
follo wing the MMSN model (Eg. [6l) . The shaded regions in Fig. [9] represent analytically estimated 
HZ (iKasting et al.lll993l . ISelsis et al.ll2007l ). Note that the positional relationship between the disk 
inner edge and the HZ is not exact because the locations of disk inner edge are uncertain. In both 
runAl and runBl, significant amount of water was delivered by planetary migration. As a result, 
the final planets are considerably "wet" except for the innermost planets that are shielded by outer 
planets. Thus, water delivery to the HZ is rather efficient around M dwarfs and the terrestrial 
planets would be rich in water. Even if the effect of type-I migration is fully retarded, water-rich 
protoplanets migrate inward by the eccentricity damping due to gravitational drag. The right panel 
of Fig. [9] shows that water-rich planets in HZs may be usually formed for relatively slow migration. 



3.5. Dependence on Disk Mass 



So far, we have adopted = 1 for protoplanetary disk in Eq. ([T]) which is relatively large for 
the disks around = 0.2Mq stars. The reduced computational cost due to the high fg allows us 
to carry out large enough number of runs for the statistical arguments. Here we discuss how the 
results can be chan ged if we consid e r less -massive disks for M dwarfs with fg ~ 0.2, which may 
be averaged values (jRavmond et al.l (j2007l ) adopted fg ~ 0.15). As described above, the results of 
A'^-body simulations are explained well by using the timescales derived in section 2.3. We discuss 
results of planetary formation in less massive disk by applying the timescales for smaller values of 

u 



-17- 



We found that the migration speed regulates final orbital configurations of the close-in ter- 
restrial planets. The migration is slower in less massive disks. According to Eqs. (jl7p and (j26p . 
both the timescales of type-I migration and the migration induced from eccentricity damping are 
inversely proportional to the disk gas scaling factor fg. Thus, in the less massive disks, the fi- 
nal planets tend to have relatively large separations in non-resonant orbits because of the slower 
migration. 

We found that the terrestrial planets around M dwarfs are generally water-ice rich by the 
relatively fast migration of icy protoplanets due to the relatively small radius (~ 0.3AU) of the ice 
line. From Eqs. (j24p and (j27p with fg = fd = 0.2, the critical masses for retention against type-I 
migration and the migration induced from eccentricity damping are given respectively by 

Me„t,n.ig - 0.17(— j (-) [-) Me, (31) 

^ 3f ^i-x3/4//,x 3/4./ .-9/20. e X -3/2. a 

The isolation mass is (Eq. |20j ) 

Comparison of these masses shows that protoplanets just outside the ice line almost grow to the 
isolation mass before starting migration. Substituting the isolation mass into Eq. ()19p . we obtain 
the accretion time as 

Similarly, substituting Eq. (j33p into Eqs. (|17p and (|26p . we obtain the migration timescales as 

,„.(5^)-(/l)-(A)-(_^)-%_, (35, 
. ,3 . :„.(!|.)-(A)-(A)-'(_i_)-(_^)-,_, ,3.) 

where tmig is the type-I migration timescale and tdamp,a is the the timescale of the e-damping 
induced migration. Both tmig and tacc are significantly smaller than disk lifetimes of ~ 10® — 10^ 
years even for /g = /d ~ 0.2. Therefore, our finding that the close-in terrestrial planets around 
M dwarfs are rather "wet" is still valid for the averaged-mass disks with fg = fd ^ 0.2 around 
~ O.2M0. 

Note, however, that the upper limit of migration timescale, tdamp,a; is comparable to the 
disk lifetimes. If the type-I migration timescale is elongated by a factor of more than 100, tmig 
is also shorter than or comparable to the disk lifetimes. Then, the transfer of water/icy materi- 
als b y migrations of protoplanets are not efficient enough. I f the migration barrier near the ice 



line (IKretke k Linll2007l . Ilda fc Linll2008bl . iKretke et al.ll2009l ) is effective also for disks around M 
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dwarfs, the transfer i s inefficient irrespective of the reduction factor of type-I migration. Since 



Raymond et alj (|2007l ) showed that the transfer of water /icy materials by scattering is also ineffi- 



cient, the planets in HZs are not always rich in water-ice in these cases. Thus, whether the planets 
in HZs around M dwarfs are habitable or not may be strongly regulated by efficiency of type-I 
migration. 



4. CONCLUSIONS AND DISCUSSION 

We have investigated accretion of terrestrial planets from planetesimals around M dwarf stars 
through a set of iV-body simulations, including the effects of disk gas. In general, accretion of 
terrestrial planets have two stages: runway/oligarchic growth and following long-term giant impacts. 
Our simulations cover all the stages from initial 5,000 planetesimals to the final planets that are 
stable for long time after possible giant impacts, fully including gravitational interactions of all the 
bodies. 

Since M dwarfs are fainter than solar- type stars, both the HZs and ice lines are located in 
the proximity of central stars. Due to the proximity, accretion of terrestrial planets have different 
features around M dwarf than around solar- type stars that are caused by the three factors: 

(a) the effective damping by disk gas due to high gas density in inner regions, 

(b) the influence from inner protoplanets that have migrated toward the disk inner edge, 

(c) the influence from outer icy protoplanets that migrate into the terrestrial planet regions. 

Regarding factor (a), it is noticed that higher disk gas density due to the proximity overwhelms 
expected smaller disk and planetary isolation masses around M dwarfs. Around M dwarfs, the disk 
mass and consequently, isolation mass of protoplanets may be generally smaller than those around 
solar-type stars. However, due to the higher disk gas density, planet-disk tidal interactions, that 
is, eccentricity (and inclination) damping and type-I migration are more efficient than those at 
~ lAU around solar type stars. Furthermore, accretion timescale in such regions is much shorter 
than disk lifetime (~ 10^ — 10''yrs). As a result, the eccentricity damping and type-I migration 
play important roles in architecture of terrestrial planets around M dwarfs. To highlight this effect, 
for A-body simulations, we adopted disks comparable to MMSN that may be relatively massive 
among disks around M dwarfs. Even in the case without type-I migration, the migration induced 
from eccentricity damping, which is 100 times slower than type-I migration, is still fast enough to 
bring protoplanets into the terrestrial planet regions for such disks. The dominance of the effects 
of disk gas was discussed with analytical arguments (section 2.3 and 3.5). 

Regarding factor (b), resonant trapping plays an important role. Because of the efficient type-I 
migration, many protoplanets migrate toward the disk inner edge and accumulate there, usually 
trapped in mutual mean-motion resonances. We set the disk inner edge to be 0.05AU, while the 



- 19 - 



HZ is around ~ 0.1 AU, which is only a factor 2 difference. In the slow migration case, in which 
resonant trapping is so efficient, trapped protoplanets line up through almost all the simulation 
regions (0.05-0.4AU), while with full migration speed predicted by the linear theory, trapping is 
not so efficient that coagulation often occurs and the final planets tend to be concentrated in inner 
regions. In the former case, about 40 planets remain in the resonances before the disk gas removal. 
After the disk gas is removed, the orbits of the planets become unstable and giant impacts occur. 
As a result, widely-spaced ('^ 20rH), non-resonant, multiple planets are formed with relatively high 
eccentricities (~ 0.05) between disk inner edge and outer region. On the other hand, in the full 
migration case, the trapped planets, the number of which is about 5, are stable even after the disk 
gas removal and closely-packed (~ 5 — lOrn), resonant planets remain in the proximity of the disk 
inner edge with low eccentricities (< 0.01). 

Therefore, we conclude that the migration speed is a key factor for final orbital configuration 
of close-in terrestrial planets around M dwarfs. The close-in planets that on-going radial velocity 
surveys have discovered may support the slow migration. The three-planet system around Gl 581 
with ~ 0.3Mq is composed of planet b (Mpsini = 16Me,a = 0.041AU), c (SM®, 0.073AU), 
and d (7.5M^, 0.25AU). They are widely-spaced (Aa ~ 21 — 47rH) non-resonant planets that are 
consistent with our slow migration model, although this system was probably formed from heavier 
disk than the MMSN disk and may need the enhancement of surface density of icy materials 
near the ice line. Around solar-type stars, the radius of the ice line is larger, so the factor (c) 
may not be applied. However, the dependence of final configuration of close-in terrestrial planets 
on the migration speed can be applied to solar-type stars. The three-planet system around a 
K dwarf HD 40307 with ~ O.77M0 consi sts of planet b (M psini = 4.2M©,a = 0.047AU), 



c (6.9Me, 0.081 AU), and d (9.2Me, 0.13AU) (jMavor et al.ll2009l ). They are also widely-spaced 



(Aa ~ 17— 20rH), non-resonant planets. These configuration is explained by setB in our calculation. 

Note that the final state could be altered by other effects which we did not address. We will 
discuss the effects of random torques exerted by strong disk turbulence due to Magneto- Rotational 
Instability in a separate paper. Two mechanisms have been proposed which lead to the excitation 
of eccentricities of the bodies and additional collisions between them: (i) during type-I migration of 
a gas giant planet, its rnean m otion resonances (mainly 2:1 resonance) sweep the terrestrial planet 



region (e.g.. lZhou et al.ll2005l ). or (b) during the dispersal of the gas disk, se cular resonances cause d 



by the gas giant planet and the disk sweep through the inner orbits (e.g., iNagasawa et al.ll2005l ). 
However, since gas giants are generally rare around M dwarfs, these mechanisms do not play an 
important role for terrestrial planets around M dwarfs. 

Factor (c) is caused by the smaller radius of the ice line around M dwarfs. Less efficient 
planetesimal accretion due to the smaller disk surface density is overwhelmed by the faster accretion 
due to smaller radii of icy regions (> 0.3AU) than those around solar-type stars. Combined with 
factor (a), migration of the icy protoplanets into the terrestrial planet regions is so efficient that 
the largest final planets have significant mass of icy components except for the innermost one that 
is shielded from impacts of icy protoplanets. In the case of enhanced surface density of ice near the 
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ice line, the largest final planets are mostly composed of ice but not rocks. It is interesting that the 
estimated bulk density for a 20Mffl transiting pla. net, GJ 436b, around a M dwarf is consistent with 
water ice (jGillon et al.ll2007l . iDeming et alj 120071 ). Although the mass and number of final planets 
are affected by this factor as well as the inner boundary conditions for type-I migration and radial 
gradient of Tig (see section 3.3), the stability of the resonantly trapped planets after the disk gas 
removal, that is, orbital separation, resonant or non-resonant configuration, eccentricities of final 
planets are determined only by migration speed. 

For the disks with fg that is several times smaller than that of our fiducial model {fg = 1), 
which may be typical disks around M dwarfs, accretion and migration timescales are still much 
shorter than disk lifetime, so the formed close-in planets are abundant in water-ice. However, if 
type-I migration speed is more than 1 00 times slower than that predicted by the linear theory o r the 
migration is trapped near the ice line ( Kretke &: Lin 2007 . Ida fc LinlboOSb . Kretke et al. 2009 ). the 
final close-in terrestrial planets would be rocky due to the inefficient water de livery by the migration , 
be cause rad i al mix ing of planetesimals is also inefficient around M dwarfs (j Raymond et al.l 120071 ) . 
As iLissauerl ()2007l ) pointed out, the ice line evolves during relatively long pre- main sequence phase 
of M dwarfs. We need to carry out detailed iV-body simulations in low-mass disks, including the 
evolution of the ice line in a future work, to clarify the details on how wet are the planets in HZs 
around M dwarfs. 



Our results around M dwarfs and calculations by iTerquem &: Papaloizoul (j2007l ) around G 
dwarfs suggest that existence of close-in relatively-large terrestrial planets are robust. However, 
our solar system does not have any close-in planet inside 0.4 AU and large fraction of extraso- 
lar planetary systems may not have the close-in super-Earths either. We will address this issue 
elsewhere. 

Our conclusion is that the migration speed determines diversity of final orbital configuration of 
close-in terrestrial planets around M dwarfs through the stability of the planets trapped in mutual 
mean-motion resonances, so we will study more detailed dependence on the migration speed as well 
as the dependence on conditions of inner disk regions in a separate paper. Around M dwarfs, these 
planets could be in HZs. Characteristics of habitable planets around M dwarfs are significantly 
affected by the details of these formation mechanisms. Although these close-in planets are well 
inside HZs for solar-type stars, their formation mechanisms may be similar. Future observations of 
close-in terrestrial planets around M dwarfs as well as solar-type stars by radial velocity surveys from 
ground and transit surveys from space will constrain the migration efficiency and the quantitative 
features of inner disk regions. 
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A. EXPRESSION FOR AERODYNAMICAL GAS DRAG FORCE 



The aerodynamic gas drag force per unit mass is (jAdachi et al.lll976l ) 



1 2 
-^CBTTTpPgAuAu, 



(Al) 



where Cd = 0.5 is the gas dra g coefficient, is the physical radius of the body and M is its mass. 



The density of disk gas pg is (jHayashi 



1981 



Pg = 2.0 X 10-V, 



'11/4 



gem 



(A2) 



'^VlAU. 

where Am is the relative velocity of the body to the disk gas. Du e to pressure gradie nt the velocity 
of disk gas fgas is smaher than Kepler velocity vk by a fraction (jAdachi et al.lll976l ) 

r \i/2 



V 



VK MAU. 
where the temperature distribution of an optically thin disk given by Eq. ([2]) is used. 



(A3) 



B. EXPRESSION FOR GRAVITATIONAL GAS DRAG FORCE 



Tanaka et al.l (j2002l ) and iTanaka &: WardI (j2004l ) derived the damping forces exerted on the 
planet, through three-dimensional linear calculation, 



-f^dampjr 
-^damp,0 
-^damp,2; 



F ■ 

^ mig,r 



mig,6 



mig,z 



}L\[!^Yf^\ni2A';.[vo-rn] + A^,Vr 



C.o / 



,2 



^.17^ 



0, 



(Bl) 

(B2) 

(B3) 
(B4) 

(B5) 
(B6) 



where -Fdamp is the specific damping force for e and i, -Fmig is the specific damping force for a and 
Q is the Keplerian angular velocity. The coefficients are given by 



Aj, 



At 



0.057 
-0.868 
-1.088 



0.176 
0.325 
-0.871. 
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C. EXPRESSION FOR POST-NEWTONIAN GRAVITY FROM THE HOST 

STAR 



The spec ific force induced by the post-Newtonian gravity from the host star, -Freh given in 
Kiddeil jlddd ) is 



rel 



G(M, + M) 



x{[(l + 3^)i'2-2(2 + ^) 



G{M, +M) 3 



2 r 



(CI) 



where c is the hght speed, v is the velocity vector, v = \v\, r = dr/dt, and = M^:M / {M^^ + M)^. 

To the lowest ord er of the eccentricity, the precession rate of the periastron is expressed as 
dMardling fc Linlbood ) 

dw 3nGM^ 

dt ac^ ' 



(C2) 



where n is the mean motion. The post-Newtonian effect is important only to the short-period 
planets. 
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Fig. 1. — Time evolution of a system on the a — e plane. The circles represent bodies and the radii 
of the circles are proportional to the physical radii of the bodies. Note that overlapping circles does 
not mean the planets actually overlap each other because the size of circles are expanded so as to 
be easy to see. The system initially consists of 5000 planetesimals. The numbers of bodies are 2970 
(IOOOTk), 2144 (IOOOOTk), 1862 (20000Tk), 1448 (50000Tk), 1107 (IOOOOOTk), 390 (SOOOOOTk), 
170 (IOOOOOOTk), and 26 (SOOOOOOTk). Tk is Keplerian time at 0.1 AU around a O.2M0 star, which 
is ~ 0.071yr. In the electronic version, bodies with M > O.OIM®, M > O.IM0, and M > are 
evpressed with blue, red, and green circles respectively. 
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Fig. 2. — Evolution of semimajor axes of massive planets. At each time for runAl, 30 most massive 
planets are plotted. Tk is Keplerian time at 0.1 AU around a O.2M0 star. The circles represent 
bodies and their radii are proportional to the physical radii of the bodies. 
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Fig. 3. — The same as Fig. [T] but the effect of type-I migration is not included. The system 
initially consists of 5000 planetesimals. The numbers of bodies are 2961 (1000 Tk), 2158 (10000 Tk), 
1864 (20000Tk), 1462 (50000Tk), 1146 (IOOOOOTk), 462 (500000Tk), 231 (IOOOOOOTr), and 83 
(5000000 Tk). 
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Fig. 4. — The same plots as Fig. [2] for runBl. In this case, the effect of type-I migration in not 
included. 
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Fig. 5. — Orbital evolution of successive simulation of runAl {left) and runBl {right) in a gas-free 
environment, neglecting small bodies. In runBl, planets suffer giant impacts, which results in losing 
of commensurabilities. 
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Fig. 6. — Orbital evolution of runCl. Left: The result of iV-body simulation with the effect of disk 
gas. Right: Successive simulation of runCl after the removal of disk gas. The effect of reversed 
torque due to inverse pressure gradient near the disk edge is included. 
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Fig. 7. — The same as Fig. [6] for runDl. The enhanced ice condensation factor (r/icc = 14) is 
adopted. 




Fig. 8. — The same as Fig. [7] for runEl. The enhanced ice condensation factor (r/ice = 14) is 
adopted, while type-I migration is not included. 
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Fig. 9. — Water mass fraction of the final planets for runAl (left) and runBl (right) after gas 
dissipation. The shaded regions express the HZ around a 0.01 Lq star. Vertical lines represent the 
location of the disk inner edge (0.05 AU). 
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Table 1. Initial Conditions and Final Results for Each Run 



run 


type-I 


reversal 


Vice 


N 


M„,ax (Me) 




M„iax,dis (Me) 


Aa(rH) 


e 


MMR 


runAl 


Yes 


No 


3 


6 


0.63 


6 


0.63 


8.3 


0.0053 


Yes 


runA2 


Yes 


No 


3 


5 


0.79 


5 


0.79 


9.5 


0.0064 


Yes 


runA3 


Yes 


No 


3 


5 


0.75 


4 


1.0 


11 


0.019 


Yes 


runA4 


Yes 


No 


3 


4 


1.1 


3 


1.1 


9.3 


0.0036 


Yes 


runBl 


No 




3 


45 


0.21 


11 


0.65 


18 


0.086 


No 


runB2 


No 




3 


38 


0.24 


17 


0.41 


16 


0.031 


No 


runB3 


No 




3 


40 


0.15 


14 


0.53 


22 


0.056 


No 


runB4 


No 




3 


36 


0.21 


14 


0.42 


19 


0.048 


No 


runCl 


Yes 


Yes 


3 


5 


1.4 


3 


1.4 


5.6 


0.0084 


Yes 


runC2 


Yes 


Yes 


3 


4 


1.6 


4 


1.6 


13 


0.011 


Yes 


runDl 


Yes 


No 


14 


3 


3.4 


3 


3.4 


8.2 


0.0083 


Yes 


runD2 


Yes 


No 


14 


3 


3.4 


3 


5.7 


11 


0.010 


Yes 


runD3 


Yes 


No 


14 


4 


3.6 


4 


3.6 


10 


0.024 


Yes 


runD4 


Yes 


No 


14 


3 


4.1 


3 


4.1 


9.3 


0.0034 


Yes 


runEl 


No 




14 


13 


0.88 


12 


3.4 


18 


0.19 


No 


runE2 


No 




14 


17 


1.0 


15 


1.4 


16 


0.073 


No 


runES 


No 




14 


12 


0.94 


17 


2.3 


16 


0.020 


No 


runE4 


No 




14 


18 


0.83 


12 


2.7 


18 


0.13 


No 



Note. — Calculation conditions and final results. The second to forth columns are calculation conditions. 
The fifth and sixth columns are results of A'^-body simulations before the removal of gas. A'' is the number of 
final planets, the masses of which are > O.OlMe, and Mmax is the mass of the largest planet. The seventh 
to eleventh columns are results of the long term evolution after the gas removal. N^is is the number of final 
planets, the masses of which are > O.OlMe, and Mmax.dis is the mass of the largest planet, Aa is the average 
orbital separations of final planets and e is the average orbital eccentricity. The final column expresses 
commensurate relationships: "Yes" means that there exists at least one commensurate relationship between 
the final planets, while "No" means that there exists no commensurate relationship. 



